Datos de tiempo de supervivencia

R para Ciencia de Datos en Salud:
Análisis Descriptivo e Inferencia Estadística

Percy Soto-Becerra M.D., M.Sc(c)

InkaStats Data Science Solutions | Medical Branch
@github/psotob91

Datos de tiempo de supervivencia

El problema de los datos de tiempo a evento

A menudo tenemos interés en describir los datos relacionados al tiempo de ocurrencia de algún evento de interés.
Los eventos de interés pueden ser muchos: muerte, cura, desarrollo de alguna enfermedad, entre abandono al tratamiento, etc.
Un problema que tenemos es que no siempre podemos observar el tiempo de ocurrencia del evento:
- Los eventos pueden tardar mucho tiempo y el estudio acaba antes de que todos tengan el evento.
- Los participantes abandonan el estudio.

Situación ideal: Todos desarrollan el evento

Idealmente, todos los participantes desarrollan el evento de interés.
- No hay censuras, los datos están completos.
Podemos estimar entonces:
- El promedio de tiempo a la muerte.
  - Otras métricas como la mediana, percentiles, etc.
- También podemos estimar la incidencia acumulada (IA) de muerte para cualquier punto del tiempo de manera directa.
  - IA a 1 mes, 2 meses, 3 meses, …, 10 meses

Todos desarrollan el evento

Media de tiempo a muerte:

\[\bar{x}_{tiempo} = \frac{10+8+10+6+10}{5} = 8.8 \text{ meses}\]

IA de muerte a 6 meses (%):

\[IA(t = 6 \text{ meses}) = 100 \times \frac{1}{5} = 20\%\]

IA de muerte a 8 meses (%):

\[IA(t = 8 \text{ meses}) = 100 \times \frac{2}{5} = 40\%\]

IA de muerte a 10 meses (%):

\[IA(t = 10 \text{ meses}) = 100 \times \frac{5}{5} = 100\%\]

No todos desarrollan el evento, pero todos completan el seguimiento

Los tiempos al evento están disponibles para algunos sujetos.
- En otros sujetos, los tiempos al evento solo dan información parcial (censurados).
Promediar estos tiempos no es válido
- Realmente no tenemos los tiempos a evento completos.
Todavía podemos hallar la IA de manera directa.
- Usamos la proporción observada de evento en cada punto de tiempo de interés.

No todos desarrollan el evento, pero todos completan el seguimiento

Media con los datos observados con tiempos incompletos:
- ¡Subestimación del verdadero valor de ella!
IA de muerte a 6 meses (%):

\[IA(t = 6 \text{ meses}) = 100 \times \frac{1}{5} = 20\%\]

IA de muerte a 8 meses (%):

\[IA(t = 8 \text{ meses}) = 100 \times \frac{2}{5} = 40\%\]

IA de muerte a 10 meses (%):
- No se puede calcular porque todo acabó en 10 meses

Situación real: No todos desarrollan el evento y hay pérdidas de seguimiento

IA no siempre puede calcularse de proporción observada en cada tiempo t.
- Hay censuras por pérdida de seguimiento.
- No todos los individuos censurados tienen misma oportunidad de ser seguidos en todo el estudio.
No todos los individuos censurados tuvieron mismo tiempo de seguimiento.
- No es justo comparar un individuo sin evento con 3 meses de seguimiento con otro de 9 meses de seguimiento.

Métodos para describir datos de supervivencia

¿Qué métodos podemos usar?

Debemos usar algún método tenga en cuenta los tiempos de seguimiento desiguales en los indviduos censurados.
- Densidad de incidencia
- Tasa de incidencia
- Incidencia y Sobrevidad Acumuladas¹
- Tablas de vida²
- Hazards
Nosotros solo nos centraremos en dos medidas que se usan tradicionalmente para describir datos de supervivencia:
- Densidad/Tasa de incidencia
- Incidencia y Sobrevida Acumulada¹
- Tabla de vida¹

Se calcula dividiendo el número de casos nuevos sobre el total de tiempo de seguimiento por persona.

\[\text{Densidad de Incidencia} = \frac{\text{Número de eventos nuevos}}{\text{Total de tiempo-persona}}\]

En el ejemplo de la derecha:

\[\text{Densidad de Incidencia} = \]

\[\frac{2 \text{ casos}}{9 +8+5+1+14 \text{ meses-persona}}\]

\[= 0.04 \text{ casos por mes-persona}\]

\[= 4 \text{ casos por 100 meses-persona}\]

Medida de la frecuencia con la que un evento de interés ocurre a lo largo de un periodo de tiempo especificado.

Algunos autores diferencian el concepto de densidad del de tasa.
Skzlo, et al., en su libro “Epidemiology beyond the basics”, hace tal distinción e indica lo siguiente:
- La densidad se calcula con los datos de individuales de los tiempos de seguimiento de cada persona (tiempo en riesgo).
- La tasa se calculan cuando no se cuentan con datos de los tiempos individuales; en cambio, se disponen de datos agregados.
  - Por tanto, se hace algunos supuestos sobre el tiempo de seguimiento de los datos agregados y se calculan estas tasas.
Densidad de incidencia y tasa de incidencia tienen las misma interpretación y supuestos.

Si tenemos datos con los tiempos de seguimiento, es relativamente fácil calcular.
Veamos los datos:

head(erc_data)

     sex age crea eGFR_ckdepi grf_cat   time5y eventd5y
1 Female  72 1.60    31.67301     G3b 3.359343        0
2 Female  84 1.98    22.49146      G4 5.000000        0
3 Female  67 1.10    51.55401     G3a 5.000000        0
4   Male  75 1.36    50.47194     G3a 3.589322        0
5 Female  71 1.24    43.41346     G3b 3.764545        0
6 Female  74 1.06    51.60754     G3a 3.657769        0
               eventd5ylab death5y
1 Alive w/o Kidney Failure       0
2 Alive w/o Kidney Failure       0
3 Alive w/o Kidney Failure       0
4 Alive w/o Kidney Failure       0
5 Alive w/o Kidney Failure       0
6 Alive w/o Kidney Failure       0

Variables de interés para cálculo de densidad de incidencia:
- time5y: tiempo de seguimiento (años)
- death5y: muerte(1) o vivo aún(0)
Cálculo de densidad de incidencia:

0.0492 caso por 1 año-persona

erc_data %>% 
  summarise(DI = sum(death5y) / sum(time5y))

          DI
1 0.04924135

4.92 casos por 100 años-persona

erc_data %>% 
  summarise(DI = sum(death5y) / sum(time5y) * 100)

        DI
1 4.924135

Muchas veces, los datos individuales no tienen los tiempos de seguimiento, sino las fechas.
- Hay que calcular los tiempos de seguimiento restando las fechas de interés.

Veamos los siguientes datos:

ebola_data_small <- readRDS("ebola_data_small.rds")
head(ebola_data_small)

# A tibble: 6 × 4
  case_id date_onset date_outcome outcome
  <chr>   <chr>      <chr>        <chr>  
1 5fe599  2014-05-13 <NA>         <NA>   
2 8689b7  2014-05-13 2014-05-18   Recover
3 11f8ea  2014-05-16 2014-05-30   Recover
4 b8812a  2014-05-18 <NA>         <NA>   
5 893f25  2014-05-21 2014-05-29   Recover
6 be99c8  2014-05-22 2014-05-24   Recover

Primero, convertir fechas de character a date:

library(lubridate)
ebola_data_small <- ebola_data_small %>% 
  mutate(date_onset = ymd(date_onset), 
         date_outcome = ymd(date_outcome)) 

head(ebola_data_small)

# A tibble: 6 × 4
  case_id date_onset date_outcome outcome
  <chr>   <date>     <date>       <chr>  
1 5fe599  2014-05-13 NA           <NA>   
2 8689b7  2014-05-13 2014-05-18   Recover
3 11f8ea  2014-05-16 2014-05-30   Recover
4 b8812a  2014-05-18 NA           <NA>   
5 893f25  2014-05-21 2014-05-29   Recover
6 be99c8  2014-05-22 2014-05-24   Recover

Luego, calculamos el tiempo de seguimiento restando fechas:

ebola_data_small <- 
  ebola_data_small %>%
  mutate(tseg_dias = as.double(date_outcome - date_onset))
head(ebola_data_small)

# A tibble: 6 × 5
  case_id date_onset date_outcome outcome tseg_dias
  <chr>   <date>     <date>       <chr>       <dbl>
1 5fe599  2014-05-13 NA           <NA>           NA
2 8689b7  2014-05-13 2014-05-18   Recover         5
3 11f8ea  2014-05-16 2014-05-30   Recover        14
4 b8812a  2014-05-18 NA           <NA>           NA
5 893f25  2014-05-21 2014-05-29   Recover         8
6 be99c8  2014-05-22 2014-05-24   Recover         2

Para contar el número de eventos, necesitamos convertir la variable outcome de character a numeric:

ebola_data_small <- 
  ebola_data_small %>% 
  mutate(outcome2 = case_when(outcome == "Recover" ~ 0, 
                              outcome == "Death" ~ 1, 
                              TRUE ~ as.numeric(NA)))
head(ebola_data_small)

# A tibble: 6 × 6
  case_id date_onset date_outcome outcome tseg_dias outcome2
  <chr>   <date>     <date>       <chr>       <dbl>    <dbl>
1 5fe599  2014-05-13 NA           <NA>           NA       NA
2 8689b7  2014-05-13 2014-05-18   Recover         5        0
3 11f8ea  2014-05-16 2014-05-30   Recover        14        0
4 b8812a  2014-05-18 NA           <NA>           NA       NA
5 893f25  2014-05-21 2014-05-29   Recover         8        0
6 be99c8  2014-05-22 2014-05-24   Recover         2        0

Cálculo de la densidad de incidencia:

Recordar que hay NA

ebola_data_small %>% 
  summarise(DI = sum(outcome2) / sum(tseg_dias))

# A tibble: 1 × 1
     DI
  <dbl>
1    NA

Si hay NA, siempre hay que considerarlo

ebola_data_small %>% 
  summarise(DI = sum(outcome2, na.rm = TRUE) / sum(tseg_dias, na.rm = TRUE) * 100)

# A tibble: 1 × 1
     DI
  <dbl>
1  4.75

4.84 casos por 100 años-persona

Con datos agregados, se suele tener el número de eventos y el tiempo de seguimiento por cada unidad de observación agregada.

Veamos los datos:

tasa_inc_data

  ncasos   region pop_mid_year
1   1245     Lima       214555
2   5234 Arequipa       341421
3   2312 Ayacucho        54545
4   4524    Tacna       123454

Calculemos la tasa de incidencia para todos los datos

tasa_inc_data %>% 
  summarise(tasa_incidencia = sum(ncasos) / sum(pop_mid_year) * 1000)

  tasa_incidencia
1        18.14094

La tasa de incidencia de la enfermedad para todas las regiones estudiadas fue de 18.14 casos por 1000 años-persona.

Calculemos la tasa de incidencia para cada región

tasa_inc_data %>% 
  mutate(tasa_incidencia = ncasos / pop_mid_year * 1000)

  ncasos   region pop_mid_year tasa_incidencia
1   1245     Lima       214555        5.802708
2   5234 Arequipa       341421       15.330047
3   2312 Ayacucho        54545       42.387020
4   4524    Tacna       123454       36.645228

La tasa de incidencia de la enfermedad para Lima fue de 5.8 casos por 1000 años-persona. En Arequipa fue de 15.3 casos por 1000 años-persona. La TI de la enfermedad en Ayacucho fue de 42.4 casos por 1000 años-persona y en Tacna fue de 36.6 casos por 1000 años-persona.

Incidencia y Supervivencia Acumulada

IA
SA
Método de KM
{survival}
{ggsurvfit}
Más info

La incidencia acumulada en un tiempo t o \(IA(t)\) o proporción de incidencia en un tiempo t, es la probabilidad de desarrollar un evento nuevo en un intervalo t de tiempo.

Si no hay pérdidas de seguimiento, se calcula directamente:

\[IA(t) = \frac{\text{Nº de eventos nuevos hasta tiempo t}}{\text{Total de individuos en riesgo al incio del estudio}}\]

Si hay pérdidas de seguimiento, debe usarse algún método que tenga en cuenta estas censuras:
- El método de Kaplan y Meier es una opción popular en Bioestadística y Epidemiología.
- Otros métodos también existen, por ejemplo, el estimador de Nelson-Aalen.

La sobrevida acumulada en tiempo t o \(SA(t)\) es la probabilidad de que un individuo sobreviva más allá de un tiempo \(t\).

Si no hay pérdidas de seguimiento, también se calcularía directamente:

\[SA(t) = \frac{\text{Nº de individuos aún en riesgo en tiempo t}}{\text{Total de individuos en riesgo al incio del estudio}}\]

Con pérdida de seguimiento, se puede usar el método de KM.
La \(SA(t)\) es el complemento de la \(IA(t)\):

\[SA(t) = 1 - IA(t)\]

También conocida como método de producto-límite, permite estimar la \(SA(t)\) y la \(IA(t)\) teniendo en cuenta que hay censuras.

No entraremos en detalles técnicos.
- Esto se ve en un curso de análisis de supervivencia.
KM permite estimar cada \(SA(t)\) e \(IA(t)\) en los puntos de tiempo \(t\) obervados.
- El reporte de cada \(SA(t)\) e \(IA(t)\) en una tabla, se llama tabla de vida
Usaremos los datos de ébola:

head(ebola_data_small, 2)

# A tibble: 2 × 6
  case_id date_onset date_outcome outcome tseg_dias outcome2
  <chr>   <date>     <date>       <chr>       <dbl>    <dbl>
1 5fe599  2014-05-13 NA           <NA>           NA       NA
2 8689b7  2014-05-13 2014-05-18   Recover         5        0

No tenemos una sola variable de interés observada, tenemos un vector de dos variables:
- Tiempo de seguimiento observado:
  - tseg_dias
- Indicador de evento y censura:
  - outcome2: 1 = muerte, 0 = censura (recuperado)
Importante que el indicador de evento/censura sea numérico con 1 y 0.

Usaremos el paquete {survival} que usa principalmente sintaxis R base.

library(survival)

Debemos definir un objeto de supervivencia con función Surv():

Surv(tiempo_seguimiento, evento_censura)

Luego, pedimos a R que estime la tabla de vida con survfit() y la guarde en un objeto:

tabla_vida <- survfit(Surv(tiempo_seguimiento, evento_censura), 
                      data = datos)

Por último, pedimos que muestre la tabla de vida:

summary(tabla_vida)

Usando los datos de ebola, creamos el objeto tabla de vida y pedimos ver sus resultados con summary:

surv_ebola <- survfit(Surv(tseg_dias, outcome2) ~ 1, 
                      data = ebola_data_small)
summary(surv_ebola)

Call: survfit(formula = Surv(tseg_dias, outcome2) ~ 1, data = ebola_data_small)

2389 observations deleted due to missingness 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    1   3499      30    0.991 0.00156        0.988        0.994
    2   3464      69    0.972 0.00281        0.966        0.977
    3   3385     149    0.929 0.00435        0.920        0.937
    4   3216     194    0.873 0.00565        0.862        0.884
    5   2989     214    0.810 0.00667        0.797        0.824
    6   2747     210    0.748 0.00740        0.734        0.763
    7   2470     179    0.694 0.00790        0.679        0.710
    8   2232     167    0.642 0.00827        0.626        0.659
    9   1992     145    0.595 0.00853        0.579        0.612
   10   1772     109    0.559 0.00870        0.542        0.576
   11   1592     119    0.517 0.00885        0.500        0.535
   12   1398      89    0.484 0.00895        0.467        0.502
   13   1214      55    0.462 0.00902        0.445        0.480
   14   1097      43    0.444 0.00908        0.427        0.462
   15    989      31    0.430 0.00913        0.413        0.448
   16    895      48    0.407 0.00923        0.389        0.426
   17    775      29    0.392 0.00931        0.374        0.411
   18    698      21    0.380 0.00938        0.362        0.399
   19    616       7    0.376 0.00941        0.358        0.395
   20    554       4    0.373 0.00944        0.355        0.392
   21    497      13    0.363 0.00957        0.345        0.383
   22    427       8    0.357 0.00969        0.338        0.376
   23    378       5    0.352 0.00979        0.333        0.372
   24    346       4    0.348 0.00989        0.329        0.368
   25    305       4    0.343 0.01002        0.324        0.363
   26    274       3    0.339 0.01014        0.320        0.360
   27    242       1    0.338 0.01019        0.319        0.359
   29    190       1    0.336 0.01029        0.317        0.357
   38     60       1    0.331 0.01155        0.309        0.354

No podemos obtener directamente estimados de \(IA(t)\), estos habría que calcularlos.

El problema con el paquete survival() es que usa R base, por lo que es difícil de manipualr.
El paquete ggsurvfit() utiliza internamente survival() pero permite trabajar con R tidy potenciando el reporte de tablas de vida.

library(ggsurvfit)

Se usa la función survfit2():

tab_vida2 <- survfit2(Surv(tseg_dias, outcome2) ~ 1, data = datos)

Para el reporte, usamos tidy_survfit, el cual remplaza a summary()

tidy_survfit(tab_vida2)

Supervivencia acumulada

surv_ebola <- survfit2(Surv(tseg_dias, outcome2) ~ 1, data = ebola_data_small)
tidy_survfit(surv_ebola)

# A tibble: 59 × 14
    time n.risk n.event n.censor cum.event cum.censor estimate std.error
   <dbl>  <dbl>   <dbl>    <dbl>     <dbl>      <dbl>    <dbl>     <dbl>
 1     0   3499       0        0         0          0    1       0      
 2     1   3499      30        5        30          5    0.991   0.00157
 3     2   3464      69       10        99         15    0.972   0.00289
 4     3   3385     149       20       248         35    0.929   0.00468
 5     4   3216     194       33       442         68    0.873   0.00647
 6     5   2989     214       28       656         96    0.810   0.00823
 7     6   2747     210       67       866        163    0.748   0.00989
 8     7   2470     179       59      1045        222    0.694   0.0114 
 9     8   2232     167       73      1212        295    0.642   0.0129 
10     9   1992     145       75      1357        370    0.595   0.0143 
# … with 49 more rows, and 6 more variables: conf.high <dbl>, conf.low <dbl>,
#   estimate_type <chr>, estimate_type_label <chr>, monotonicity_type <chr>,
#   conf.level <dbl>

Incidencia acumulada

tidy_survfit(surv_ebola, type = "risk")

# A tibble: 59 × 14
    time n.risk n.event n.censor cum.event cum.censor estimate std.error
   <dbl>  <dbl>   <dbl>    <dbl>     <dbl>      <dbl>    <dbl>     <dbl>
 1     0   3499       0        0         0          0  0         0      
 2     1   3499      30        5        30          5  0.00857   0.00157
 3     2   3464      69       10        99         15  0.0283    0.00289
 4     3   3385     149       20       248         35  0.0711    0.00468
 5     4   3216     194       33       442         68  0.127     0.00647
 6     5   2989     214       28       656         96  0.190     0.00823
 7     6   2747     210       67       866        163  0.252     0.00989
 8     7   2470     179       59      1045        222  0.306     0.0114 
 9     8   2232     167       73      1212        295  0.358     0.0129 
10     9   1992     145       75      1357        370  0.405     0.0143 
# … with 49 more rows, and 6 more variables: conf.high <dbl>, conf.low <dbl>,
#   estimate_type <chr>, estimate_type_label <chr>, monotonicity_type <chr>,
#   conf.level <dbl>

Tiempos específicos de interés

tidy_survfit(surv_ebola, times = c(1, 2, 5), type = "risk")

# A tibble: 3 × 14
   time n.risk n.event n.censor cum.event cum.censor estimate std.error
  <dbl>  <dbl>   <dbl>    <dbl>     <dbl>      <dbl>    <dbl>     <dbl>
1     1   3499      30        5        30          5  0.00857   0.00157
2     2   3464      69       10        99         15  0.0283    0.00289
3     5   2989     557       81       656         96  0.190     0.00823
# … with 6 more variables: conf.high <dbl>, conf.low <dbl>,
#   estimate_type <chr>, estimate_type_label <chr>, monotonicity_type <chr>,
#   conf.level <dbl>

Selecciona solo columnas de interés

tidy_survfit(surv_ebola,  times = c(1, 2, 5), type = "risk") %>% 
  select(time, n.event, n.risk, estimate)

# A tibble: 3 × 4
   time n.event n.risk estimate
  <dbl>   <dbl>  <dbl>    <dbl>
1     1      30   3499  0.00857
2     2      69   3464  0.0283 
3     5     557   2989  0.190

Más información sobre análisis de datos de supervivencia:

{survival}: https://cran.r-project.org/web/packages/survival/index.html (Ver Sección Vignettes)
{ggsurvfit}: https://www.danieldsjoberg.com/ggsurvfit/index.html
Un tutorial básico para epidemiólogos sobre análisis de supevivencia: https://epirhandbook.com/en/survival-analysis.html
Una lista amplia de paquetes en R para análisis de supervivencia: https://cran.r-project.org/web/views/Survival.html
Un (más avanzado) tutorial interesante sobre análisis de supervivencia: https://www.emilyzabor.com/tutorials/survival_analysis_in_r_tutorial.html

Curvas de supervivencia

Definición
{survival}
{ggsurvfit}
{survminer}
Más info

Es la versión gráfica de la tabla de vida.

Hay varios métodos para obtenerla.
- No paramétricos: método de KM, Nelson-Aalen, etc.
- Semiparamétricos: predichas de modelo de regresión de Cox.
- Paramétricos: predichas de modelos paramétricos de supervivencia.
Veremos algunas funciones para obtener gráficos para reporte.
Solo nos centraremos en el método no paramétrico de Kaplan-Meier para describir datos de supervivencia.

Se puede graficar el objeto de survfit con la función plot.

surv_ebola <- survfit(Surv(tseg_dias, outcome2) ~ 1, 
                      data = ebola_data_small)
summary(surv_ebola)

Call: survfit(formula = Surv(tseg_dias, outcome2) ~ 1, data = ebola_data_small)

2389 observations deleted due to missingness 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    1   3499      30    0.991 0.00156        0.988        0.994
    2   3464      69    0.972 0.00281        0.966        0.977
    3   3385     149    0.929 0.00435        0.920        0.937
    4   3216     194    0.873 0.00565        0.862        0.884
    5   2989     214    0.810 0.00667        0.797        0.824
    6   2747     210    0.748 0.00740        0.734        0.763
    7   2470     179    0.694 0.00790        0.679        0.710
    8   2232     167    0.642 0.00827        0.626        0.659
    9   1992     145    0.595 0.00853        0.579        0.612
   10   1772     109    0.559 0.00870        0.542        0.576
   11   1592     119    0.517 0.00885        0.500        0.535
   12   1398      89    0.484 0.00895        0.467        0.502
   13   1214      55    0.462 0.00902        0.445        0.480
   14   1097      43    0.444 0.00908        0.427        0.462
   15    989      31    0.430 0.00913        0.413        0.448
   16    895      48    0.407 0.00923        0.389        0.426
   17    775      29    0.392 0.00931        0.374        0.411
   18    698      21    0.380 0.00938        0.362        0.399
   19    616       7    0.376 0.00941        0.358        0.395
   20    554       4    0.373 0.00944        0.355        0.392
   21    497      13    0.363 0.00957        0.345        0.383
   22    427       8    0.357 0.00969        0.338        0.376
   23    378       5    0.352 0.00979        0.333        0.372
   24    346       4    0.348 0.00989        0.329        0.368
   25    305       4    0.343 0.01002        0.324        0.363
   26    274       3    0.339 0.01014        0.320        0.360
   27    242       1    0.338 0.01019        0.319        0.359
   29    190       1    0.336 0.01029        0.317        0.357
   38     60       1    0.331 0.01155        0.309        0.354

plot(surv_ebola)

Este gráfico no es {ggplto2}, en cambio, es de R base.
- Se puede personalizar, pero la sintaxis de estos gráficos es relativamente complicada y más limitada que la sintaxis de {ggplot2}

Podemos usar la función ggsurvfit() para realizar gráficos publicables de la KM.
- Estos gráficos sí siguen sintaxis {ggplot2} y uno puede seguir agregando capas.

La función ggsurvfit() es un wrapper de ggplot() con geom_step() diseñado para crear la curva de supervivencia:

surv_ebola %>% 
  ggsurvfit()

Le podemos agregar una tabla de riesgo, la cual en Bioestadística es un estándar mínimo de publicación de curvas de KM:

surv_ebola %>% 
  ggsurvfit() + 
  add_risktable()

También se puede agregar marcas de censuras:

surv_ebola %>% 
  ggsurvfit() + 
  add_risktable() + 
  add_censor_mark()

Puedes configurar temas

surv_ebola %>% 
  ggsurvfit() + 
  add_risktable() + 
  add_censor_mark() + 
  theme_classic()

También temas de la tabla de riesgos

surv_ebola %>% 
  ggsurvfit() + 
  add_risktable(theme = theme_risktable_boxed()) + 
  add_censor_mark() +
  theme_classic()

Puedes hacer más personalizaciones:

surv_ebola %>% 
  ggsurvfit(size = 1, color = "red") + 
  add_risktable(theme = theme_risktable_boxed()) + 
  add_censor_mark(color = "red") + 
  theme_classic() + 
  labs(y = "Supervivencia Acumulada (%)", 
       x = "Tiempo (Días)", 
       title = "Curva de supervivencia en pacientes con Ébola")

El paquete {survminer} también ofrece funciones para gráficos ggplot2().
- Sin embargo, los objetos que crean también tienen atributos adicionales:
  - No usa la lóigica de capas ni el símbolo +
- Personalmente creo que esto genera algunas inconsistencias y obliga a salirse un poco de la sintaxis {ggplot2}
No veremos {survminer} en detalle, porque ggsurvfit() lo supera por creces.
Solo veremos un par de ejemplos para conocer, por cultura general, el paquete en general:

library(survminer)

El gráfico de curvas KM puede obtener mediante:

surv_ebola %>% 
  ggsurvplot()

Puedes agregar una tabla de riesgo :

surv_ebola %>% 
  ggsurvplot(risk.table = TRUE)

También puedes hacer algunas personalizaciones si son argumentos dentro de la función:

surv_ebola %>% 
  ggsurvplot(risk.table = TRUE, 
             xlab = "Tiempo (días)", 
             ylab = "Supervivencia Acumulada (%)", 
             ggtheme = theme_minimal(), 
             risk.table.title = "Número en riesgo")

Más información sobre gráficos KM para datos de supervivencia:

{ggsurvfit}: https://www.danieldsjoberg.com/ggsurvfit/index.html
- Gráficos: https://www.danieldsjoberg.com/ggsurvfit/articles/gallery.html
- Temas: https://www.danieldsjoberg.com/ggsurvfit/articles/themes.html
- Listado de todas las funciones: https://www.danieldsjoberg.com/ggsurvfit/reference/index.html
{survminer}: https://rpkgs.datanovia.com/survminer/index.html
Un tutorial para epidemiólogos sober survminer: https://epirhandbook.com/en/survival-analysis.html

Nuestro turno

Abra el proyecto var_tiempo.Rproj y dentro de este, abra el archivo quarto var_tiempo.qmd.
Siga las instrucciones indicadas en este.
Renderice el archivo quarto final.

10:00